﻿using System;
using System.Collections.Generic;
using System.Text;

namespace ParallelMath
{
    public static class NonParallelMath
    {

        public static void Multiply(float[,] R, float[,] A, float[,] B)
        {
            int
              aRows = A.GetLength(0),
              aCols = A.GetLength(1),
              bRows = B.GetLength(0),
              bCols = B.GetLength(1),
              rRows = R.GetLength(0),
              rCols = R.GetLength(1);
            if (aCols == bRows && rRows == aRows && rCols == bCols)
            {
                for (int r = 0; r < rRows; r++)
                    for (int c = 0; c < rCols; c++)
                    {
                        float sum = 0.0f;
                        for (int k = 0; k < aCols; k++)
                            sum += A[r, k] * B[k, c];
                        R[r, c] = sum;
                    }
            }
        }

        public static double ludcmp(float[,] a, int n, int[] index)
        {
            int i, imax = 0, j, k;
            float d, big, dum, sum, temp;
            float[] vv = new float[n];

            d = 1.0f;
            for (i = 0; i < n; i++)
            {
                big = 0.0f;
                for (j = 0; j < n; j++)
                {
                    if ((temp = Math.Abs(a[i, j])) > big)
                    {
                        big = temp;
                    }
                }
                if (big == 0.0)
                {
                    
                    throw new Exception("Singular matrix in routine ludcmp");
                }
                vv[i] = 1.0f / big;
            }

            for (j = 0; j < n; j++)
            {
                for (i = 0; i < j; i++)
                {
                    sum = a[i, j];
                    for (k = 0; k < i; k++)
                    {
                        sum -= a[i, k] * a[k, j];
                    }
                    a[i, j] = sum;
                }
                big = 0.0f;
                for (i = j; i < n; i++)
                {
                    sum = a[i, j];
                    for (k = 0; k < j; k++)
                    {
                        sum -= a[i, k] * a[k, j];
                    }
                    a[i, j] = sum;
                    if ((dum = vv[i] * Math.Abs(sum)) >= big)
                    {
                        big = dum;
                        imax = i;
                    }
                }
                if (j != imax)
                {
                    for (k = 0; k < n; k++)
                    {
                        dum = a[imax, k];
                        a[imax, k] = a[j, k];
                        a[j, k] = dum;
                    }
                    d = -d;
                    vv[imax] = vv[j];
                }
                index[j] = imax;
                if (a[j, j] == 0.0)
                {
                    a[j, j] = 1.0e-20f;
                }
                if (j != n - 1)
                {
                    dum = 1.0f / a[j, j];
                    for (i = j + 1; i < n; i++)
                    {
                        a[i, j] *= dum;
                    }
                }
            }

            return d;
        }




        public static float[,] inverse(float[,] a)
        {
            int i, j, N = a.GetUpperBound(0) + 1;
            int[] index = new int[N];
            float[] col = new float[N];

            float[,] y = new float[N, N];

            ludcmp(a, N, index);

            for (j = 0; j < N; j++)
            {
                for (i = 0; i < N; i++) col[i] = 0.0f;
                col[j] = 1.0f;
                lubksb(a, N, index, col);
                for (i = 0; i < N; i++)
                {
                    y[i, j] = col[i];
                }
            }

            return y;
        }

        public static void lubksb(float[,] a, int n, int[] index, float[] b)
        {
            int i, ii = -1, ip, j;
            float sum;

            for (i = 0; i < n; i++)
            {
                ip = index[i];
                sum = b[ip];
                b[ip] = b[i];
                if (ii != -1)
                {
                    for (j = ii; j <= i - 1; j++)
                    {
                        sum -= a[i, j] * b[j];
                    }
                }
                else if (sum != 0.0)
                {
                    ii = i;
                }
                b[i] = sum;
            }

            for (i = n - 1; i >= 0; i--)
            {
                sum = b[i];
                for (j = i + 1; j < n; j++)
                {
                    sum -= a[i, j] * b[j];
                }
                b[i] = sum / a[i, i];
            }
        }
        
    }
}

    

